Using R as a Geographic Information System

# Query CMS Provider Data Catalog

query <- paste0(
  "https://data.cms.gov/provider-data/api/1/datastore/query/nrdb-3fcy/0",
  "?conditions[0][property]=State&conditions[0][value]=", params$state,
  "&conditions[0][operator]=="
)
  
raw <- jsonlite::fromJSON(query)
  
hospital_df <- raw$results

# Identify hospitals w/out L&D

no_inpatient <- hospital_df |>  
  filter(stringr::str_detect(score, "does not provide inpatient labor")) |> 
  distinct(facility_id, .keep_all = TRUE) |> 
  select(c(facility_id, facility_name)) 

# Preview 

head(hospital_df)

head(no_inpatient)

# Filter and reformat hospital data

hospital_df <- hospital_df |> 
  distinct(facility_id, .keep_all = TRUE) |> 
  mutate(inpatient = case_when(facility_id %in% no_inpatient$facility_id  ~ "No",
                               .default = "Yes")) |> 
  mutate(full_address = paste0(address, ", ", citytown, ", ", state, " ", zip_code),
         .after = zip_code) |> 
  select(-c(telephone_number:end_date))

# Geocode hospital addresses

hospital_df <- hospital_df |> 
  tidygeocoder::geocode(address = full_address, method = "arcgis") 

# Convert to sf object

hospitals <- hospital_df |>
  st_as_sf(coords = c("long", "lat"), 
           crs = 4326) # use {leaflet} preferred CRS

# Preview for quick QC

mapview::mapview(hospitals)

# Set HERE API key

set_key(Sys.getenv("HERE_API_KEY"))

# Loop over point data to make isochrones file (using Hannah Recht's function)

## select hospitals w/ L&D

ip_hospitals <- hospitals |> 
  filter(inpatient == "Yes")

## create empty vectors for loop output

isochrones <- vector(mode = "list", length = nrow(ip_hospitals))
error_rows <- vector(mode = "list", length = nrow(ip_hospitals))

for (i in 1:nrow(ip_hospitals)) {
  
    print(i)
  
    ## get isochrone for that point, using delay to avoid rate limiting
    Sys.sleep(0.9)
    
    ## filter to ith point
    point_temp <- ip_hospitals %>% filter(row_number() == i)
    point_id <- point_temp$facility_id
    
    isochrones_temp <- tryLocation(point_temp)
    
    ## save any errored out points
    if (is.null(isochrones_temp)) {
        error_rows <- bind_rows(error_rows, point_temp)
    } else {
        isochrones <- bind_rows(isochrones, isochrones_temp)    
    }
}

# Remove extraneous columns

isochrones <- isochrones |> 
  select(facility_id, geometry)

# Left join hospital data

ip_hospitals <- ip_hospitals |> 
  select(c(facility_id, facility_name, full_address, inpatient)) |> 
  st_drop_geometry()

isochrones <- isochrones |> 
  left_join(ip_hospitals, by = join_by(facility_id))

# Reproject

isochrones <- isochrones |> 
  st_transform(crs = 4326)

# Preview

mapview::mapview(isochrones)

# Download county boundaries 

counties <- counties(state = c(params$state), 
                     cb = TRUE, 
                     resolution = "500k", 
                     year = 2023) 

# Reproject to {leaflet} preferred CRS

counties <- counties |> 
  st_transform(crs = 4326)

# Create a quick QC viz

mapview::mapview(counties)

# Specify text for popups

hosp_popup <- paste0(
  "<b>", hospitals$facility_name, "</b>", "<br/>",
  hospitals$full_address
)

iso_popup <- paste0(
  "This is a 45-minute driving radius around ", isochrones$facility_name, "."
)

# Create color palette

color_icon <- awesomeIcons(
  icon = 'location-dot',
  iconColor = 'black',
  library = 'fa',
  markerColor = ifelse(hospitals$inpatient == 'Yes', 'green', 'gray')) # no HEX codes

# Make an interactive map

leaflet() %>%
  addProviderTiles("CartoDB.Voyager") %>%    # add base map
  addPolygons(data = counties,
              weight = 1,
              color = "#999999",
              fillColor = "#d4d4d4",
              stroke = TRUE,
              popup = paste0(counties$NAME)) %>%
  addAwesomeMarkers(data = hospitals, 
             icon = color_icon,
             group = hospitals$facility_name,
             popup = ~hosp_popup) %>%
  addPolygons(data = isochrones,
              weight = 1, opacity = 1.0,     # set stroke width and opacity
              stroke = TRUE,
              color = "#669933",
              fillColor = "#669933", 
              fillOpacity = 0.5, 
              group = isochrones$facility_name,
              popup = ~iso_popup,
              highlightOptions = highlightOptions(color = "#527c29",
                                                  weight = 3,
                                                  bringToFront = TRUE)
  ) %>%
  addLayersControl(
    overlayGroups = sort(c(isochrones$facility_name, hospitals$facility_name)),
    options = layersControlOptions(collapsed = TRUE)
  )